function params = CalcMapParams(dir)

% a set of parameters of direction (orientation) selectivity are obtained by vector averaging.
% dir: input stack of direction or orientation dF images (one image per condition)
% th: preferred angle obtained by vector averaging, normalized to 0-1.
% mag: vector magnitude
% ave_change: average signal change to all directions
% max_change: max signal change to any directions
% tune: vector sum / scalar sum of singal changes to all directions, according to Bonhoeffer et al. (1995)
% this function can be used to estimate parameters of both orientation & direction selectivity.
% if input is dF, mag, ave_change and max_change will be absolute change.
% if input is ratio(dF/F), mag, ave_change and max_change will be ratio change.
% th & tune will be the same, regardless the input is dF or ratio (dF/F).
% when the signal change to one direction is negative, it is replaced by zero, because tune becomes strange.
% thus, tune parameter is between 0 and 1.


dim = size(dir);
xsize = dim(1);
ysize = dim(2);
ndir = dim(3);

params.th = zeros(xsize, ysize);
params.mag = zeros(xsize, ysize);
params.ave_change = zeros(xsize, ysize);
params.max_change = zeros(xsize, ysize);
params.tune = zeros(xsize, ysize);

a = zeros(ndir,1);

for x = 1:xsize
    for y = 1:ysize
        Vx = 0;
        Vy = 0;
        sum = 0;
        % vector averageing
        for i = 1:ndir
            if dir(x,y,i) < 0; 
                a(i) = 0;
            else
                a(i) = dir(x,y,i); 
            end;
            Vx = Vx + a(i) * cos(2 * (i - 1) * pi / ndir);
            Vy = Vy + a(i) * sin(2 * (i - 1) * pi / ndir);
            sum = sum + a(i);
        end
        params.th(x,y) = atan2(Vy, Vx);
        params.mag(x,y) = (Vx ^ 2 + Vy ^2 ) ^ 0.5;
        params.ave_change(x,y) = sum ./ ndir;
        params.max_change(x,y) = max(a);
        if sum~=0; 
            params.tune(x,y) = params.mag(x,y) / sum; 
        else
            params.tune(x,y) = 0;
        end;
    end
end

% params.th = params.th./pi/2+0.5; % normalize to 0-1.
params.th = mod(params.th, 2 * pi)/2/pi; % the other line is not consistent with \analysis\VectorAverage.m, which does the mod "flip" of the negative values to the upper half of the range